module StochasticSimulation.GillespieAlgorithm

open System
open StochasticSimulation.ChemicalReaction.ChemicalReaction
open StochasticSimulation.ChemicalReaction.ChemicalSpecies
open System.Collections.Generic

(* Gillespie Algorithm

	1. initialise
	2. calcaulte propensity a_i for each reaction, total propensity a0
	3. generate randoms r1 and r2
	4. calcaulte timestep = (1/a0)*ln(1/r1)
	5. Determine which reaction has occurred
	5. update counts
	6. t <- t + timestep
*)

    
let chooseReaction propensities total r2  =
    let points = propensities |> Seq.scan (fun rt amount -> (rt + snd amount) ) 0.0 |> Seq.skip 1
    let partitions = Seq.zip (Seq.map (fun x -> x/total) points ) propensities    
    Seq.pick (fun (x,p) -> if r2 < x then Some(fst p) else None) partitions
     
let executeNextReaction (reactions:seq<Reaction>) (t:float) (rnd:Random) (system:seq<SpeciesCount>) =
    let (r1,r2) = (rnd.NextDouble(),rnd.NextDouble())
    let propensities = reactions |> Seq.map ( fun r -> (r, r.Propensity system) )
    let a0 = propensities |> Seq.fold (fun acc (r,a) -> acc + a ) 0.0
    let timestep = (log (1.0/r1)) / a0
    let nextReaction = chooseReaction propensities a0 r2
    let updatedSystem = nextReaction.React(system)
    { tau =(t + timestep) ; UpdatedSystem = Seq.toList updatedSystem }